@yuzhesong A few more comments:

NS::PulsarAccretion() needs @return in description
NS::PulsarAccretion() parameter kappa should be p_Kappa
NS::PulsarAccretion() has a couple of stray std::cout ... statements - should be removed

NS::UpdateMagneticFieldAndSpin() has a bunch of stray std::cout ... statements - should be removed
NS::UpdateMagneticFieldAndSpin() indentation is out somewhere - lines 635 and 636 have braces in the same column
NS::UpdateMagneticFieldAndSpin() what is the comment at the end of line 599 for? There is no P_f in NS::UpdateMagneticFieldAndSpin()...
NS::UpdateMagneticFieldAndSpin() ditto line 625
NS::UpdateMagneticFieldAndSpin() the comment at the end of line 619 is not useful - remove it

We use camelCase in COMPAS for variable and function names, not snake_case, and we (generally) start local variable names with lowercase, and function names with uppercase. In NS::UpdateMagneticFieldAndSpin(), change:
Accretion_Results to accretionResults
Next_Accretion_Results to nextAccretionResults
divide_timestep_by to divideTimestepBy
do_not_stop_looping to doNotStopLooping
new_timestep_size to newTimestepSize
new_massGain to newMassGain

In NS::UpdateMagneticFieldAndSpin(), this block of code doesn't need to be inside the for loop:

if (n==1) {
    last_B = std::get<0>(Accretion_Results);
    last_f = std::get<1>(Accretion_Results);
    last_fdot = std::get<2>(Accretion_Results); 
    last_am = std::get<3>(Accretion_Results);
}

What is the motivation for setting Next_Accretion_Results = Accretion_Results ; at line 572, when you break out of the for loop at line 573, and Next_Accretion_Results is not used outside that for loop?

Can I suggest you rewrite the existing while loop something like this:

int intraSteps = 100;
bool done      = false;
while (!done) {
    
    intraSteps *= 2;
    
    double thisTimestepSize = p_Stepsize / intraSteps;
    double thisMassGain     = p_MassGainPerTimeStep * 1000.0 / intraSteps; 
    
    accretionResults = PulsarAccretion(m_PulsarDetails.magneticField, m_PulsarDetails.spinFrequency, 
    m_AngularMomentum_CGS, thisTimestepSize, thisMassGain, p_Kappa, p_Epsilon);

    B    = std::get<0>(accretionResults);
    f    = std::get<1>(accretionResults);
    fdot = std::get<2>(accretionResults); 
    am   = std::get<3>(accretionResults);

    int count = 0;
    while (!done) {
        thisAccretionResults = PulsarAccretion(B, f, am, thisTimestepSize, thisMassGain, p_Kappa, p_Epsilon);
        f  = std::get<1>(thisAccretionResults);
        
        if (utils::Compare(f, 0.0) < 0) break;

        B    = std::get<0>(thisAccretionResults); 
        fdot = std::get<2>(thisAccretionResults); 
        am   = std::get<3>(thisAccretionResults); 
    
        if (++count >= intraSteps) done = true;        
    }
}

I'll stop here for now and come back to it later.

EDIT:

I forgot to ask - is fdot used anywhere in the code (or going to be)? If not, you don't need to extract it from accretionResults.

And can we put in a guard for the while loop never terminating? Maybe we do the inner loop (or the outer loop - doesn't really matter I guess) a maximum of X times (you pick X)?

You can also use std::tie() to receive the results from a function that returns a tuple:

std::tie(B, f, fdot, am) = PulsarAccretion(B, f, am, thisTimestepSize, thisMassGain, p_Kappa, p_Epsilon);

@jeffriley the last two commits should've addressed all of your comments.

@SimonStevenson and I have been doing some testing on the new MT descriptions (#1125 ) by @ilyamandel and check its impact on neutron star in mass transfer. The figure attached here plots, from top to bottom, using the configuration file attached:

• P-Pdot diagram of a NS in a NS-MS binary system
• Period of NS vs. time
• Pdot of NS vs. time
• Mass of NS vs. time
• Mdot of NS vs. time,
• magnetic field of NS vs. time.

The blue curves (labeled I) in all figure are produced with MT treatments in #1125 , and orange curves (labeled I+S) are produced with implementation in #1125 AND Simon’s treatment on breaking down the nuclear mass transfer of a MS onto NS into smaller timesteps, which is currently implemented in this PR.

Treatments in #1125 still has a one short step of large mass transfer, while the other treatments provide gradual mass transfer. In the end, both treatments transfer the same amount of mass onto the NS, but within different timescale (ass seen in the plot). This means that magnetic field of the NS should decay to the same value (as seen in the plot). But due to the different Mdot (as shown in the second last panel of the plot), the Alfven radius is then different, which then cause the difference in the calculation of P and Pdot of the pulsar at each time step, giving different final results. Overall, it seems like having the mass transfer broken down over time seems a better treatment for pulsar evolution.

Simon and I think that maybe we can use Ilya’s mass transfer treatments in all cases, unless we have evolve-pulsars=True, then we can implement Simon's treatment as a special case for NS in MT?

compasConfig_v2.txt

Hi @yuzhesong ,
Unfortunately, #1125 was not the complete fix I hoped, due to some underlying inconsistencies elsewhere in the code. I hope that #1132 will address the remaining issue.

Just in case, to help me with testing your particular issue, could you give me a couple of COMPAS command lines (i.e., initial conditions for the binaries at ZAMS) that should lead to slow MT onto NSs?

@ilyamandel will wait for #1132 to be merged to check again.

If you use the compacConfig_v2.txt file listed in my reply above, and change it to a yaml file, it should produce a NS-MS binary system in slow case A mass transfer.

The mass of the NS grows from 1.1773 Msol at formation to 1.4322 Msol at the end of the evolution in your example with the latest PR... but then the binary merges, which should not happen. I'll move the PR to draft and investigate ASAP, but probably not before Monday or even Tuesday.

@ilyamandel that is weird and shouldn't have happened. Using the corrections in #1125 and Simon's treatments, this binary system evolves to "allowed time exceeded".

Actually, managed to fix that faster than I thought. You'll be happy to know that we now do, in fact, limit nuclear timescale MT to no more than the MT rate times the timestep duration. However, we only do this for nuclear timescale mass transfer, and we do it for all accretors, so I think this is a more robust solution. I have tested that your NS accretes around ~0.2 Msol (slowly, over many timesteps). I have not tested other features such as the B field strength -- leave that for you, after PR #1132 is reviewed.

@jeffriley @ilyamandel @SimonStevenson

Ilya's new mass transfer fixes work.
Updated the codes for this PR for review to merge to dev.

@yuzhesong

"Updated the codes for this PR for review to merge to dev." I don't see any updates - am I missing them? Or are they older changes?

Also, there are conflicts to resolve. And finally - the PR is marked as a draft - do you want to mark it ready for review?

@jeffriley Sorry I didn't manage to push the changes made to PR properly. It should all be resolved now hopefully.

@ilyamandel working on cleaning up the code.

For your question regarding line 601:
double Jacc = MoI * PPOW(G_CGS * mass_g /r_cm_3, 0.5) * p_MassGainPerTimeStep * 1000.0 / mass_g ;

We are calculating the added angular momentum here as

J_acc = dJ = (J / M) * dm = MoI * sqrt(GM/R^3) * dm/M

====================
Also regarding your suggestion of using an off-the-shelf ODE integrator: I'm not particularly familiar with this tool, so I will have to look into implementing one for this part of the code. Do you have any other resources regarding this topic?

Hi @yuzhesong :

General info on numerical ODE integration that uses higher-order integrators to improve accuracy/cost ratios and automatically chooses step sizes to achieve desired tolerances: Numerical Recipes textbook.

Specifically, I used Boost ODE integrators; there's documentation online (see https://www.boost.org/doc/libs/1_82_0/libs/numeric/odeint/doc/html/index.html ), but you are probably best off just using my example in BaseBinaryStar::CalculateMassTransferOrbit().

J_acc = dJ = (J / M) * dm = MoI * sqrt(GM/R^3) * dm/M

Are you assuming that the star is rotating at the Keplerian orbital frequency prior to accretion (that's what J = MoI * sqrt(GM/R^3) implies)?

Closing at @yuzhesong 's suggestion since this is replaced by #1315 .